The woylier package implements tour interpolation paths between frames using Givens rotations. This provides an alternative to the geodesic interpolation between planes currently available in the tourr package. Tours are used to visualise high-dimensional data and models, to detect clustering, anomalies and non-linear relationships. Frame-to-frame interpolation can be useful for projection pursuit guided tours when the index is not rotationally invariant. It also provides a way to specifically reach a given target frame. We demonstrate the method for exploring non-linear relationships between currency cross-rates.
When data has up to three variables, visualization is relatively intuitive, while with more than three variables, we face the challenge of visualizing high dimensions on 2D displays. This issue was tackled by the grand tour (Asimov 1985) which can be used to view data in more than three dimensions using linear projections. It is based on the idea of rotations of a lower-dimensional projection in high-dimensional space. The grand tour allows users to see dynamic low-dimensional (typically 2D) projections of higher dimensional space. Originally, Asimov’s grand tour presents the viewer with an automatic movie of projections with no user control. Since then new work has added interactivity to the tour, giving more control to users (Buja et al. 2005). New variations include the manual (Cook and Buja 1997) or radial tour (Laa et al. 2023), little tour, guided tour (Cook et al. 1995), local tour, and planned tour. These are different ways of selecting the sequence of projection bases for the tour, for an overview see Lee et al. (2022).
The guided tour combines projection pursuit with the grand tour and it is implemented in the tourr package (Wickham et al. 2011). Projection pursuit is a procedure used to locate the projection of high-to-low dimensional space that should expose the most interesting feature of data, originally proposed in Kruskal (1969). It involves defining a criterion of interest, a numerical objective function that indicates the interestingness of each projection, and an optimization for selecting planes with increasing values of the function. In the literature, a number of such criteria have been developed based on clustering, spread, and outliers.
A tour path is a sequence of projections and we use an interpolation to produce small steps simulating a smooth movement. The current implementation of tour in tourr package uses geodesic interpolation between planes. The geodesic interpolation path is the locally shortest path between planes with no within-plane spin (see Buja et al. (2004) for more details). As a result, the rendered target plane could be a within-plane rotation of the target plane originally specified. This is not a problem when the structure we are looking for can be identified from any rotation. However, even simple associations in 2D, such as the calculated correlation between variables, can be very different when the basis is rotated.
Most projection pursuit indexes, particularly those provided by the tourr are rotationally invariant. However, there are some projection pursuit index where the orientation of frames does matter. One example is the splines index proposed by Grimm (2016). The splines index computes a spline model for the two variables in a projection, in order to measure non-linear association. It can be useful to detect non-linear relationships in high-dimensional data. However, its value will change substantially if the projection is rotated within the plane (Laa and Cook 2020). The procedure in Grimm (2016) was less affected by the orientation because it considered only pairs of variables, and it selects the maximum value found when exchanging which variable is considered as predictor and response variable.
Figure 1 illustrates the rotational invariance problem for a modified spline2D index, where we always consider the horizontal direction as the predictor variable, and the vertical direction as the response. Thus, our modified index computes the splines on one orientation, exaggerating the rotational variability. The example data was simulated to follow a sine curve and the modified splines index is calculated on different within-plane rotations of the data. Although they have the same structure, the index values vary greatly.
The lack of rotation invariance of the splines index raises complications in the optimisation process in the projection-pursuit guided tour as available in the tourr. Fixing this is the motivation of this work. The goal with the frame-to-frame interpolation is that optimisation would find the best within-plane rotation, and thus appropriately optimize the index.
Figure 1: The impact of rotation on a spline index that is NOT rotation invariant. The index value for different within-plane rotations take very different values: (a) original projection has maximum index value of 1.00, (b) axes rotated 45\(^o\) drops index value to 0.83, (c) axes rotated 60\(^o\) drops index to a very low 0.26. Geodesic interpolation between planes will have difficulty finding the maximum of an index like this because it is focused only on the projection plane, not the frame defining the plane.
A few alternatives to geodesic interpolation were proposed by Buja et al. (2005) including the decomposition of orthogonal matrices, Givens decomposition, and Householder decomposition. The purpose of the woylier package is to implement the Givens paths method in R. This algorithm adapts the Given’s matrix decomposition technique which allows the interpolation to be between frames rather than planes.
This article is structured as follows. The next section provides the theoretical framework of the Givens interpolation method followed by a section about the implementation in R. The method is applied to search for nonlinear associations between currency cross-rates.
The tour method of visualization is animated high-to-low dimensional data rotation that is a movie, one-parameter (time) family of static projections. Algorithms for such dynamic projections Buja et al. (2005) are based on the idea of smoothly interpolating a discrete sequence of projections.
The topic of this article is the construction of paths of projections. Interpolation of paths of projection can be compared to connecting line segments that interpolate points in Euclidean space. Interpolation acts as a bridge between continuous animation and discrete choice of sequences of projections.
The interpolating paths of plane versus frames
Current implementation of tourr package is locally shortest (geodesic) interpolation of planes. The pitfall of this interpolation method is that it does not account for rotation variability. Therefore, the interpolation of frames is required when the orientation of projection matters. If the rendering on a frame and on the rotated version of the frame yields the same visual scenes, it means the orientation does not matter.
The orientation of frames could be important when non-linear projection pursuit function is used in guided tour. An illustration of such cases are shown in Figure 2.
Figure 2: Plane to plane interpolation (left) and Frame to frame interpolation (right). We used dog index for illustration purposes. For some non-linear index orientation of data could affect the index.
Before continuing with the interpolation algorithms, we need to define the notations.
Let the \(p\) be the dimension of original data and \(d\) be the dimension onto which the data is being projected.
A frame \(F\) is defined as \(p\times d\) matrix with pairwise orthogonal columns of unit length that satisfies, where \(I_d\) is the identity matrix in d dimensions.
\[F^TF = I_d\]
Paths of projections are given by continuous one-parameter families \(F(t)\) where \(t\in [a, z]\) interval representing time. We denote the starting frame by \(F_a = F(a)\) and target frame by \(F_z = F(z)\). Usually, \(F_z\) is selected target basis that has chosen via various methods. While grand tour chooses target frames randomly, guided tour chooses the target plane by optimizing the projection pursuit index. Interpolation methods are used to move from \(F_a\) to \(F_z\).
\(B\) is preprojection basis of \(F_a\) and \(F_z\).
Preprojection algorithm
In order to make the interpolation algorithm simple, we need to carry out “preprojection” step. The purpose of preprojection is to limit data subspace that the interpolation path, \(F(t)\), is traversing. In other words, preprojection step make sure the interpolation path between two frames \(F_a\) and \(F_z\) is not going to the data space that is not related to \(F_a\) and \(F_z\). Simply, prepojection algorithm is defining the joint subspace of \(F_a\) and \(F_z\).
The procedure starts with forming an orthonormal basis by applying Gram-Schmidt to \(F_z\) with regard to \(F_a\). We denote this orthonormal basis by \(F_\star\). Then build preprojection basis \(B\) by combining \(F_a\) and \(F_\star\) as follows:
\[B = (F_a, F_{\star})\]
The dimension of the resulting orthonormal basis, \(B\), is \(p\times 2d\).
Then, we can express the original frames in terms of this basis:
\[F_a = B^TW_a, F_z = B^TW_z\]
The interpolation problem is then reduced to the construction of paths of frames \(W(t)\) that interpolate the preprojected frames \(W_a\) and \(W_z\). Because \(B\) is orthonormalized basis of \(F_z\) with regard to \(F_a\), \(W_a\) is \(2d\times d\) matrix of 1, 0s. This is an important character for our interpolation algorithm of choice, Givens interpolation.
Givens interpolation path algorithm
A rotation matrix is a transformation matrix used to perform a rotation in Euclidean space in a plane. A rotation matrix that transforms 2D plane by an angle \(\theta\) looks like this:
\[ \begin{bmatrix}\cos \theta &-\sin \theta \\\sin \theta &\cos \theta \end{bmatrix} \]
If the rotation is in the plane of selected 2 variables, it is called Givens rotation. Let’s denote those 2 variables \(i\) and \(j\). The Givens rotation is useful for introducing zeros on a grand scale and used for computing the QR decomposition of matrix in linear algebra problems. One advantage over other transformation methods which is particularly useful in our case is the ability to zero elements more selectively.
The interpolation methods in the woylier package is based on the fact that in any vector of a matrix, one can zero out the \(i\)-th coordinate with a Givens rotation (Golub and Loan 1989) in the \((i, j)\)-plane for any \(j\neq i\). This rotation affects only coordinate \(i\) and \(j\) and leave all other coordinates unchanged.
Sequences of Givens rotations can map any orthonormal d-frame F in p-space to standard d-frame \(E_d=((1, 0, 0, ...)^T, (0, 1, 0, ...)^T, ...)\).
The interpolation path construction algorithm from starting frame \(F_a\) to target frame \(F_z\) is illustrated below. The example is 2D path construction process of original 6D data frame.
In our example, \(F_a\) and \(F_z\) are \(p\times d\) or \(6\times2\) matrices that are orthonormal. The preprojection basis \(B\) is \(p\times 2d\) matrix that is \(6\times 4\).
In our example, \(W_a\) looks like:
\[ \begin{bmatrix}1 & 0 \\0 &1 \\ 0&0 \\0&0\end{bmatrix} \]
\(W_z\) is orthonormal \(2d\times d\) matrix that looks like:
\[ \begin{bmatrix} a_{11} & a_{12} \\a_{21} &a_{22} \\ a_{31}&a_{32} \\a_{41}&a_{42}\end{bmatrix} \]
\[ W_a = R_m(\theta_m) ... R_2(\theta_2)R_1(\theta_1)W_z\]
At each rotation, the angle \(\theta_i\) that zero out the second coordinate of a plane is calculated.
When \(d = 2\), there are 5 rotations involved with 5 different angles that makes each elements 0. For example, the first rotation angle \(\theta_1\) is an angle in radian between \((1, 0)\) and \((a_{11}, a_{21})\). This rotation matrix would make element \(a_{21}\) zero:
\[R_1(\theta_1) = G(1, 2, \theta_1) = \begin{bmatrix} cos\theta_1 & -sin\theta_1 & 0 & 0 \\sin\theta_1 &cos\theta_1 & 0 &0 \\ 0&0&1&0 \\0&0&0&1\end{bmatrix}\]
6th rotation is not necessary due to orthonormality of columns. If we make one element of a column 1 that means all other elements must be 0.
\[R(\theta) = R_1(-\theta_1) ... R_m(-\theta_m), \ W_z = R(\theta)W_a\]
Performing these rotations would go from the starting frame to the target frame in one step. But we want to do it sequentially in a number of steps so interpolation between frames looks dynamic.
Next step should include the time parameter, \(t\), so that it shows the interpolation process rendered in the movie-like sequence. We break \(\theta_i\) into the number of steps, \(n-step\), that we want to go from starting frame to the target frame, which means it moves by equal angle in each step.
Finally, we reconstruct our original frames using \(B\). This reconstruction is done at each step of interpolation so that we have interpolated path as result. We use \(F_t\) to project the orignal data into lower dimensions.
\[F_t = B W_t\]
Projection pursuit index functions
The properties of several projection pursuit index functions were
investigated in Laa and Cook (2020). The smoothness,
squintability, flexibility, rotation invariance, and speed of projection
pursuit index functions were examined. The one property that is
interesting to us is rotation invariance. The rotational invariance is
examined by computing projection pursuit index for different rotations
within 2D plane. It is established that the dcor2d,
splines2d and TIC index are not strictly
rotationally invariant. The splines2d index measures
nonlinear association between variable by fitting spline model. It
compares the variance of residuals and the functional dependence is
stronger when the index value is larger.
We implemented each steps in Givens interpolation path algorithm in
separate functions and combined them in givens_full_path()
function. Here is the input and output of each functions and it’s
descriptions.
| name | description | input | output |
|---|---|---|---|
givens_full_path(Fa, Fz, nsteps)
|
Construct full interpolated frames. | Starting and target frame (Fa, Fz) and number of steps | An array with nsteps matrix. Each matrix is interpolated frame in between starting and target frames. |
preprojection(Fa, Fz)
|
Build a d-dimensional pre-projection space by orthonormalizing Fz with regard to Fa. | Starting and target frame (Fa, Fz) | B pre-projection p x 2D matrix |
construct_preframe(Fa, B)
|
Construct preprojected frames. | A frame and the pre-projection p x 2D matrix | Pre-projected frame in pre-projection space |
row_rot(a, i, k, theta)
|
Performs Givens rotation . | A frame and the pre-projection p x 2D matrix | theta angle rotated matrix a |
calculate_angles(Wa, Wz)
|
Calculate angles of required rotations to map Wz to Wa. | Preprojected frames (Wa, Wz) | Names list of angles |
construct_moving_frame(Wt, B)
|
Reconstruct interpolated frames using pre-projection. | Pre-projection matrix B, Each frame of givens path | A frame of on a step of interpolation |
The interface of tour is that it renders one projection of data at a time. It displays one projection and asks for the next projection. Therefore, path of projections shown below is sequence of projections to be renders at tour display.
The givens_full_path() function returns the intermediate
interpolation step projections in given number of steps. The code chunk
below demonstrates the interpolation between 2 random basis in 5
steps.
set.seed(2022)
p <- 6
base1 <- tourr::basis_random(p, d=2)
base2 <- tourr::basis_random(p, d=2)
base1
[,1] [,2]
[1,] 0.24406482 -0.57724655
[2,] -0.31814139 0.06085804
[3,] -0.24334450 0.38323969
[4,] -0.39166263 0.01182949
[5,] -0.08975114 0.59899558
[6,] -0.78647758 -0.39657839
base2
[,1] [,2]
[1,] -0.64550501 -0.17034478
[2,] 0.06108262 0.87051018
[3,] -0.03470326 0.26771612
[4,] -0.05281183 0.25452167
[5,] -0.43004248 0.27472455
[6,] -0.62502981 0.03560765
givens_full_path(base1, base2, nsteps = 5)
, , 1
[,1] [,2]
[1,] 0.02498501 -0.57102411
[2,] -0.26080833 0.26278410
[3,] -0.19820064 0.40434178
[4,] -0.35542927 0.08341593
[5,] -0.14433023 0.57626698
[6,] -0.86308174 -0.31951242
, , 2
[,1] [,2]
[1,] -0.1909937 -0.5290164
[2,] -0.1874044 0.4550600
[3,] -0.1459678 0.4046873
[4,] -0.2970111 0.1522888
[5,] -0.2003186 0.5261305
[6,] -0.8824688 -0.2197674
, , 3
[,1] [,2]
[1,] -0.38527579 -0.4457635
[2,] -0.10411664 0.6258684
[3,] -0.09577045 0.3811614
[4,] -0.22183655 0.2089977
[5,] -0.26412984 0.4533801
[6,] -0.84414115 -0.1137724
, , 4
[,1] [,2]
[1,] -0.54115467 -0.32350096
[2,] -0.01855341 0.76518462
[3,] -0.05630484 0.33422743
[4,] -0.13748432 0.24504604
[5,] -0.34020920 0.36617868
[6,] -0.75431619 -0.02150119
, , 5
[,1] [,2]
[1,] -0.64550501 -0.17305000
[2,] 0.06108262 0.86649508
[3,] -0.03470326 0.26851774
[4,] -0.05281183 0.25487107
[5,] -0.43004248 0.27511042
[6,] -0.62502981 0.03766958
In this section, we illustrate the use of givens_full_path() function by plotting the interpolated path between 2 frames. This also a way of checking if interpolated path is moving in equal size at each step.
For plotting the interpolated path of projections, we used geozoo package (Schloerke 2016). 1D projection is plotted on unit sphere, while 2D projection is visualized on torus. The points on the surface of sphere and torus shape are randomly generated by functions from the geozoo package.
Interpolated paths of 1D projection
1D projection of data in high dimension linear combination of data that is normalized. Therefore, we can plot the point on the surface of a hypersphere. Figure 3 shows the Givens interpolation steps between 2 points, 1D projection of 6D data that is.